Week 04 - Network Travel and Accessibility

Why compute travel times on a network?

  • Because you can ….
  • Because it matters….

Humans are animals whose activities are constrained by a spatial range. Some of us like to live close to work, some to shopping, others to hotspots of social activity or green areas. While people differ in specific priorities, they all tend to optimize the cost of movement. In this tutorial, based on Kyle Walker’s MUSA class and mapboxAPI package, we learn to use mapboxAPI to calculate ease of travel across modern road networks so that you see how your environment constrains and impacts your daily behavior.

Libraries

To get started with mapboxapi, you’ll need to first install some packages. mapboxapi was released to CRAN in 2020, so we can install with install.packages(). Beware that the process can be a bit entailed and require several updates. That is at least what it took me on a W10 machine running R 4.1.1 The result is worth it, but allow for enough time and restart if R hangs up in the installation/update process. It is also good to install or update the Rcpp library to >1.0.7 version.

# If running this the first time, uncomment these lines

# install.packages('Rcpp') install.packages('fasterize')
# install.packages('mapboxapi', dependencies = TRUE)

MapboxAPI

Before we get started using Mapbox services in R, you’ll need a valid Mapbox account with an access token. Fortunately, Mapbox has generous allowance of 100,000 free API requests so you should not need to pay anything for the API use. To set up your account, visit https://account.mapbox.com/auth/signup/ to establish an account - all you need to provide is an email address to sign up! Fill out the form and verify your account through the email Mapbox sends you; you’ll be taken directly to your Mapbox account dashboard page.

Note the “default public token” that appears on your dashboard screen - you’ll come back to this page in a moment. First, look to the right side of your screen and click “View billing” under Plan. This is where Mapbox will handle your billing information. Nothing you’ll do today will be intensive enough to incur charges - but if you plan to do serious work, you need to be aware that mapboxapi is not a forever free service. Copy the access token that appears on your screen to your clipboard, then return to R.

All features in mapboxapi require a valid Mapbox access token to work. Now that you have yours in hand, you can set yours up! Load the mapboxapi package and install your token as follows:

my_token <- "YOUR TOKEN GOES HERE"

library(mapboxapi)
mb_access_token(my_token, install = TRUE)

The optional argument install = TRUE saves the token to your .Renviron, allowing you to use mapboxapi functions in the future without having to worry about setting your token. To use this feature, restart your R session.

Mapbox maps

As you might expect, one of the features of Mapbox services is its ability to create neat web maps. Why do we care when we can play with leaflet() and tmap() and other libraries? It is important as the Mapbox Terms of Service require that Mapbox API outputs be visualized on Mapbox maps and we will be creating MapboxAPI outputs, such as isochrones, accessibility and navigation tools.

Mapbox maps are accessed through styles, which are custom design configurations applied to OpenStreetMap or even user-generated vector map tilesets. You’ll learn how to create and use your own map style with Mapbox later in this workshop. However, Mapbox provides a number of their styles to all users with a Mapbox access token. The most recent versions of these styles (as of the workshop date) are as follows:

  • streets-v11: The core Mapbox Streets basemap
  • outdoors-v11: A basemap designed for outdoor recreation uses
  • light-v10: A light, greyscale background suitable for thematic overlay
  • dark-v10: A dark basemap suitable for thematic overlay
  • satellite-v9: A global satellite basemap derived from MODIS, Landsat, & proprietary imagery sources
  • satellite-streets-v11: The satellite basemap with a streets overlay

You already know one of the most popular R packages for interactive data visualization in R, the Leaflet package maintained by RStudio, which wraps the Leaflet JavaScript library for web mapping. mapboxapi provides a convenience function, addMapboxTiles(), to help you use Mapbox styles in Leaflet in a straightforward way. This function queries the Mapbox Static Tiles API and converts a Mapbox style into static tiles for web mapping.

Let’s load the leaflet and mapboxapi libraries and set up an interactive map: Please note: your map won’t show up in the RStudio Viewer pane; pop it out to a web browser to view it. It should show in the rmarkdown in its own time

library(leaflet)
library(mapboxapi)

mapbox_map <- leaflet() %>% addMapboxTiles(style_id = "streets-v11", username = "mapbox")

mapbox_map

Geocoding with MapboxAPI - Single locations

Once we’ve set up our Leaflet map with a Mapbox basemap, we’ll likely want to focus it on a specific location. mapboxapi includes functionality for R users to interact with the Mapbox Search API. Implemented functions include mb_geocode() for forward geocoding, which refers to the conversion of a description of a place (like an address) into longitude/latitude coordinates; and mb_reverse_geocode(), which converts coordinates into a place description.

Both functions default to using the mapbox.places API endpoint, which is to be used for temporary geocoding. This means that the endpoint cannot be used to store geocoded information nor can it be used for batch geocoding (e.g., a spreadsheet of addresses). These tasks are permissible with the mapbox.places-permanent endpoint, which is not included with free accounts. R users looking for free batch-geocoding solutions should use other packages like the tidygeocoder or opencage package. Mapbox geocoding with the mapbox.places endpoint can be used to focus web maps and guide navigation services, which will be illustrated in the following sections.

Let’s use mb_geocode() to identify the coordinates representing the Aarhus University Nobelparken campus.

nobel <- mb_geocode("Nobelparken, Jens Chr.Skous Vej 5, Aarhus, 8000, Denmark ")
nobel
[1] 10.20558 56.17236

By default, mb_geocode() returns a length-2 vector representing the longitude and latitude coordinates of the geocoded location. The function can also return an sf POINT object or an R list representing the full API response, if requested. Using the returned coordinates, we can focus our Leaflet Mapbox map with the setView() function:

mapbox_map %>% setView(lng = nobel[1], lat = nobel[2], zoom = 14)

Quick Exercise

Try it out! Make a Leaflet map in R using a Mapbox basemap of your choice, focused on a location of your choice. For locations in non-English-speaking countries: mb_geocode() has a language argument that can be used to improve the accuracy of queries in languages other than English. Supported languages (and how to specify them) are found in the Mapbox documentation here.

## A QUICK MAP OF YOURS

Drawing isochrones with Mapbox and R

Creating and visualizing isochrones is beautifully straightforward with the mb_isochrone() function in mapboxapi. Supported travel profiles include driving (with no traffic), cycling, and walking. mb_isochrone() by default returns a simple features polygon object that can be used for visualization and even spatial analysis.

Let’s try drawing isochrones around the Nobelparken. mb_isochrone() accepts an an input a coordinate pair, a location description as a character string, or an sf object. We can use our nobel object here to initialize the isochrones around campus.

nobel_iso_drive <- mb_isochrone(nobel, profile = "driving", time = c(4, 8, 12))

nobel_iso_drive
Simple feature collection with 3 features and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 10.09455 ymin: 56.12393 xmax: 10.29138 ymax: 56.25838
geographic CRS: WGS 84
# A tibble: 3 x 2
   time                                                                 geometry
  <int>                                                            <POLYGON [°]>
1    12 ((10.12558 56.25838, 10.12547 56.25346, 10.12427 56.25336, 10.12403 56.~
2     8 ((10.17258 56.21238, 10.17197 56.21236, 10.17255 56.21233, 10.17243 56.~
3     4 ((10.21858 56.18912, 10.21799 56.18594, 10.21698 56.18575, 10.21664 56.~

An sf object is returned with a time column representing the travel-time around the location. time is organized in descending order to ensure that overlapping isochrones are plotted correctly, with the shortest time visualized last (on top).

Using Leaflet’s addPolygons() function, we can add the isochrones to our map

colors <- viridisLite::viridis(3)

# You can set your own colors, too library(RColorBrewer) colors <-
# c('#fc8d59','#ffffbf','#91bfdb')


leaflet() %>% addMapboxTiles(style_id = "outdoors-v11", username = "mapbox") %>% 
    addPolygons(data = nobel_iso_drive, color = rev(colors), fillColor = rev(colors), 
        fillOpacity = 0.5, opacity = 0.5, weight = 0.2) %>% addLegend(labels = c(4, 
    8, 12), colors = colors, title = "Drive-time<br/>around Nobel")

### Routing with mapboxapi

mapboxapi can also be used to quickly represent and visualize routes between two locations, or alternatively along multiple locations. At its simplest, however, mb_directions() just requires an origin and a destination:

route <- mb_directions(origin = nobel, destination = "Moesgaard Museum, Aarhus, Denmark", 
    profile = "cycling")

mapbox_map %>% addPolylines(data = route, popup = paste0("Distance (km): ", round(route$distance, 
    1), "<br/>Time (minutes): ", round(route$duration, 1)))

The optional argument steps = TRUE will break the route object into separate rows for each leg of the trip, and return travel instructions in a number of different languages (English is the default).

route_dir <- mb_directions(origin = nobel, destination = "Moesgaard Museum, Aarhus, Denmark", 
    profile = "cycling", steps = TRUE)

route_dir
Simple feature collection with 37 features and 3 fields
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: 10.20415 ymin: 56.08815 xmax: 10.2258 ymax: 56.17236
geographic CRS: WGS 84
First 10 features:
                         geometry distance   duration
1  LINESTRING (10.20558 56.172...   0.0195 0.22833333
2  LINESTRING (10.20527 56.172...   0.0159 0.24500000
3  LINESTRING (10.20526 56.172...   0.0462 0.85833333
4  LINESTRING (10.20492 56.171...   0.0490 0.27500000
5  LINESTRING (10.20416 56.171...   0.0037 0.03666667
6  LINESTRING (10.20415 56.171...   0.1644 0.91500000
7  LINESTRING (10.20534 56.170...   0.3655 1.46500000
8  LINESTRING (10.20736 56.167...   0.1916 0.76833333
9  LINESTRING (10.20836 56.165...   0.4981 2.32166667
10 LINESTRING (10.21144 56.161...   0.1069 0.42833333
                             instruction
1  Head west on Jens Christian Skous Vej
2             Turn left onto Nobelparken
3            Keep right onto Nobelparken
4                             Turn right
5                              Turn left
6                     Make a slight left
7     Go straight onto Nørrebrogade (O1)
8                    Make a slight right
9     Go straight onto Nørrebrogade (O1)
10    Make a slight right onto Nørreport

Exercise 1 - Navigation

Now that you’ve learned how to use isochrone and routing services in mapboxapi, try them out for yourselves! Create the following maps:

  1. An isochrone map around a location of your choice. Times can be specified at 1-minute intervals all the way up to 60 minutes using a vector.
  2. A route between two locations of your choice, using a travel profile of your choice.
## YOUR CODE HERE

How many and which Viking monuments are reachable within each of the three bands?

## YOUR CODE HERE

Analyzing Accessibility - Skejby Hospital

Resource accessibility is a major issue in modern society. While equality of access is a fundamental tenet of most democracies, the realities of geography, differential biology and income collude to prevent such equitable access. Democracies concerned for their citizens thus need to interrogate service availability to everyone regardless of location, social status and other social markers and compensate for its deficiencies.

Our class is taking place only a few months after the 2021 communal and regional elections in Denmark, where a number of issues touched on accessibility of medical care. One of the major issues on the ballot was the centralisation of hospitals in the Midtjylland region. The argument for centralisation being the better quality and greater efficiency of services for inpatients, the con being the excessive travel times for citizens residing in more remote parts of Jutland.

Hospital accessibility can be analyzed using Mapbox services and the mapboxapi package in sf-based data science workflows.

In this exercise we can explore one or more advanced applications of mapboxapi. We’ll probably not get beyond 2 but let’s see :)

  1. We’ll examine how to visualize accessibility to a Skejby hospital in Aarhus;
  2. Identify areas where populations may have difficulty reaching medical services in central Jutland; and
  3. Build a routing app with Shiny that identifies the closest hospital to a user’s address. This section may include some new concepts or techniques - but it is designed to illustrate where you can go with mapboxapi in your work!

Visualizing the (in)accessibility of Skejby

We can visualize this situation in central Jutland with layered isochrones. We already used this technique to show multiple drive times around the Aarhus University earlier in this tutorial. In this case, we will use mb_isochrone() to generate dozens of isochrones, then visualize them simultaneously to illustrate an accessibility gradient in the region.

We’ll first generate the isochrones using a vector of times, 1 through 45 at 1-minute intervals, around the Skejby hospital emergency room address.

library(mapboxapi)

isos <- mb_isochrone(location = " Palle Juul-Jensens Blvd. 161, 8200 Aarhus", profile = "driving", 
    time = 1:45)

Next, we can visualize our overlapping isochrones. We’ll use the viridis color palette as we did previously in the tutorial, and generate a color palette derived from the time column in our dataset. Once specified, we can add these polygons to our Mapbox basemap with a mostly-transparent fill opacity.

pal <- colorNumeric("viridis", isos$time, na.color = "transparent")

mapbox_map %>% addPolygons(data = isos, fillColor = ~pal(time), stroke = FALSE, fillOpacity = 0.1) %>% 
    addLegend(values = isos$time, pal = pal, title = "Drive-time to Skejby hospital")

The result illustrates some of the wide differences in accessibility between various parts of the region. One notable issue with this visualization approach, however, is that the layering of isochrones makes it difficult to view the basemap beneath them. This can be resolved by converting to a raster dataset and generating an “accessibility surface” for improved visualization.

Let’s save the isochrone layers

To make continuation easier and save on the API calls.

# dir.create('data') library(sf) write_sf(isos, 'data/isos4326.shp')

Making an “accessibility surface”

Accessibility surfaces are commonly used in geographic information systems applications to identify the distance from any particular location to a geographic feature of interest. We can apply this concept to network-based accessibility by using mapboxapi tools. To create the accessibility surface, we will convert our isochrones to a raster dataset using the fasterize package. Raster datasets represent geographic information as grid cells defined by a cell size. Higher-resolution raster datasets are represented with smaller cell sizes.

To generate the accessibility surface raster, we will need to apply a coordinate system transformation to “project” our data to two-dimensional coordinates. This will allow us to specify the raster’s resolution in meters. We generate a 100m resolution raster, and use the fasterize() function to allocate the minimum overlapping value from our isochrones to each grid cell. The result can then be mapped with Leaflet’s addRasterImage() function.

library(Rcpp)
library(fasterize)
library(raster)

isos_proj <- st_transform(isos, 25832)

template <- raster(isos_proj, resolution = 100)
iso_surface <- fasterize(isos_proj, template, field = "time", fun = "min")

pal <- colorNumeric("viridis", isos$time, na.color = "transparent")

mapbox_map %>% addRasterImage(iso_surface, colors = pal, opacity = 0.5) %>% addLegend(values = isos$time, 
    pal = pal, title = "Drive-time to Skejby Hospital")

And again, let’s save the accessibility surface for later

# library(raster) ?writeRaster() writeRaster(iso_surface100,
# 'data/iso_surface100.tif', format = 'GTiff' )

Ok, task one down. Onto task two.

Analyzing Accessibility - Hospitals in Midtjylland

So the accessibility of Skejby looks pretty good, plus it is inside a big city which has a lot of public transport options. What about more rural regions of Midtjylland and their access to medical care?

You can find a list of all hospitals in Denmark in https://www.listsclub.com/hospitals-in-denmark/. For Midtjylland, the list is:

  • Regionshospitalet Brædstrup in Brædstrup
  • Regionshospitalet Grenaa in Grenaa
  • Regionshospitalet Hammel Neurocenter in Hammel
  • Regionshospitalet Herning in Herning
  • Regionshospitalet Holstebro in Holstebro
  • Regionshospitalet Horsens in Horsens
  • Regionshospitalet Kjellerup in Kjellerup
  • Regionshospitalet Lemvig in Lemvig
  • Regionshospitalet Odder in Odder
  • Regionshospitalet Randers in Randers
  • Regionshospitalet Ringkøbing in Ringkøbing
  • Regionshospitalet Samsø on the island of Samsø
  • Regionshospitalet Silkeborg in Silkeborg
  • Regionshospitalet Skanderborg Sundhedscenter in Skanderborg
  • Regionshospitalet Skive in Skive
  • Regionshospitalet Tarm in Tarm
  • Regionshospitalet Viborg in Viborg
  • Aarhus Universitetshospital Skejby in Aarhus

Use the opencage library to geocode these. Opencage is an API service, so you will again need to sign up in order to get an API Key for 2500 free API requests per day at 1 per second rate. Installation and registration of opencage is much less entailed than mapboxapi and should work straitforwardly.

Opencage API

You have read about geocoding above in the mapboxapi section, and learnt that mapboxapi does not let you batch-geocode addresses. Opencage does with familiar-sounding oc_forward_df() and oc_backward_df() which take a single address or a vector of string addresses! The next steps will take you through the opencage API package to geocode a list of hospitals in order to get their coordinates.

  • Register at opencagedata.com/users/sign_up.
  • Generate an API key at the OpenCage dashboard.
  • Save your API key as an environment variable like OPENCAGE_KEY=yourkey in .Renviron. See help(oc_config) for alternative ways to set your OpenCage API key.

If you want to know more, the opencage vignettes are really great! There are “Introduction to opencage” vignette(“opencage”); “Customise your query” vignette(“customise_query”) and “Output options” vignette(“output_options”)

library(opencage)
library(tidyverse)

# see configuration options
help(oc_config)

# set key in the environment
Sys.setenv(OPENCAGE_KEY = "YOUR KEY")

# set the key interactively if it is missing.
oc_config()

Geocoding hospitals

Now you are ready to turn place names into latitude and longitude coordinates:

b <- oc_forward_df(placename = "Regionshospitalet Brædstrup in Brædstrup", countrycode = "DK")
leaflet() %>% addTiles() %>% addMarkers(b$oc_lng, b$oc_lat)

Let’s try with a vector of addresses. You might need to clean the list above in regex101.com before you make it into a character list! Alternatively, check out this

hospitals <- c("Regionshospitalet Brædstrup in Brædstrup", "Regionshospitalet Grenaa in Grenaa", 
    "Regionshospitalet Hammel Neurocenter in Hammel", "Regionshospitalet Herning in Herning", 
    "Regionshospitalet Holstebro in Holstebro", "Regionshospitalet Horsens in Horsens", 
    "Regionshospitalet Kjellerup in Kjellerup", "Regionshospitalet Lemvig in Lemvig", 
    "Regionshospitalet Odder in Odder", "Regionshospitalet Randers in Randers", "Regionshospitalet Ringkøbing in Ringkøbing", 
    "Regionshospitalet Samsø on the island of Samsø", "Regionshospitalet Silkeborg in Silkeborg", 
    "Regionshospitalet Skanderborg Sundhedscenter in Skanderborg", "Regionshospitalet Skive in Skive", 
    "Regionshospitalet Tarm in Tarm", "Regionshospitalet Viborg in Viborg", "Aarhus Universitetshospital Skejby in Aarhus")

Batch-geocode!

h <- oc_forward_df(placename = hospitals, countrycode = "DK", limit = 1, no_annotations = TRUE)

# Look at the object
h

# Look at the map
leaflet() %>% addTiles() %>% addMarkers(h$oc_lng, h$oc_lat, popup = h$oc_formatted)

# Save object write_csv(h, 'data/midty_hospitals.csv')

Well, the list is not perfect in all instances. Sometimes it gets the hospital , other times it places a marker in the middle of a town. If you want to fix it, provide street addresses and re-geocode.

Let’s convert the tibble to an sf POINT object. This dataset can be used to analyze which areas are immediately covered by hospital services, and which are not. We’ll measure accessibility using isochrones as above, and consider a 20 minute walk-time, cycling- and driving-time around each hospital. mb_isochrone() can accept sf objects as input, and will retain an ID from the input sf object if the column name is specified.

# Convert to sf object
h_sf <- st_as_sf(h, coords = c("oc_lng", "oc_lat"), crs = 4326)  #%>% 
# st_transform(crs = 25832)


# Walking range
library(opencage)
walking_isos <- mb_isochrone(h_sf, profile = "walking", time = 20, id = "placename")

# Cycling range
cycling_isos <- mb_isochrone(h_sf, profile = "cycling", time = 20, id = "placename")

# Driving range
driving_isos <- mb_isochrone(h_sf, profile = "driving", time = 20, id = "placename")

Accessibility of Midtjylland hospitals

These results can be visualized on our Mapbox map:

mapbox_map %>% addPolygons(data = walking_isos, popup = ~id) %>% addPolygons(data = cycling_isos, 
    color = "red") %>% addPolygons(data = driving_isos, color = "yellow")

Well done! Now you can see the reachable area within a 20-minute walk, 20-minute bike ride and a 20-min drive. You probably don’t want to bike for 20 minutes to reach a hospital if not feeling well. Likewise, car driving is comfortable, but not every place is within 20-mins drive of a hospital and not everybody has a car. Finally, remember that these isochrones do not include traffic updates.

Questions:

1. Which areas are not covered by 20-min drive or 20-min bikeride? How big is their total area?

2. How much do you need to increase the driving range so that everyone in Midtjylland can access a hospital (let’s assume everyone has a car and just want to inform them of how long it will take)?

Assessing Risk with Demographic data

There are also individuals with disabilities, the elderly, or others requiring accessible transport, who do not have access to a car. Getting to these hospital sites may prove difficult in areas outside the isochrones. Within large cities, this problem may be solved with public transport, such as the letbanen in Aarhus. Accessibility may also be less of an issue in areas where car ownership is widespread. We can analyze this additional variable with demographic data, also obtained within R.

Wrangle data from DS

In the next exercise, you will find out which kommunes face the problem of having lots of high-risk households, where there is a no car, and they are outside the isochrones and in rural areas (out of the reach of public transport)

Number of households without cars are available via Danmark Statistik via the “Familiernes bilrådighed (faktiske tal) efter rådighedsmønster, område” dataset. The “familier uden biler” sits in csv format as a no_cars.csv with the numbers corresponding to 2020 and 2021 respectively. Overall, there were 1178935 households without a car in 2020 and 1166668 in 2021. The number of course drops with time, while number of car-owning households go up steadily from 1,896,387 in 2020 to 1,929,511 in 2021

  • Load the dataset
  • rename the V3 and V4 column to nc2020 and nc2021. (In the provided dataset the numeric columns are number of carless households in year 2020 and in 2021
no_cars <- _________("../data/no_cars.csv", header = FALSE)
names___________ <- c("municipality","nc2020", "nc2021")

It would be nice to know the percentage of households without a car Next, get municipalities of Midtjylland and join up with the no_cars attributes, selecting which kommunes have the least cars and are therefore most at risk.

Percentages of inhabitants without a car per municipality

Create a simple feature or spatial polygons dataset that contains information on the percentage of households that do not have access to a car.

  • First read in the type and count of households (“husstande-og-familier”) from 2020 and 2021 from DS . Beware that the dataset has no headers, and read.csv() might perform better than tidyverse counterpar.
  • Summarize the hh object after grouping by municipality to get a sum of households per kommune. Rename the resulting columns to municipality and households and label the object hh_mun
  • Join the hh_mun to no_cars object by the shared municipality column with inner_join(). Look up the function to know how exactly to specify the shared columns.
  • In the resulting hh_nocar object, create two new columns where you calculate the percentage of carless househols for 2020 and 2021, using the total household count and the nc2020 and nc2021 columns.
  • Load population data by municipality into pop object in order to compare and contrast with households. How many people on average per household?

– You can either use the provided file (population by municipality from the first quarter of 2022) called munic_pop_2022Q1.csv or, – Concoct your own dataset (you can further break it down by sex or age) from https://www.statbank.dk/FOLK1A

  • Load municipality data either via getData() or from a saved object from Week03
  • Create no_cars_pct_mun object by joining the households with no car hh_nocar object to the municipalities with inner_join(). Beware to preserve the spatial dimension of the munic for later use in a map. Look up the function to know how exactly to specify the shared columns. You may also want to rename the pop columns first so they are easier to map to munic
## YOUR CODE HERE
[1] 5456264
[1] 11746840
  • Create a thematic map of no_cars_pct_mun by filtering municipalities that have 25+ percent of households without a car. I recommend mapview(), but you are free to use your favorite library. Mapview usually works well at the end of tidyverse pipeline but as it may interfere with other libraries, I recommend saving the script and upon restart, running only the mapview() code chunk so as to avoid interference.
  • Which areas in Midtjylland have 25 and more percent of carless households?

Now of course, you have the whole situation in Denmark. In the next steps, you need to constrain the municipalities with carless households spatially to those in Midtjylland and check the overlap with hospital isochrones, taking the area outside these isochrones with st_difference().

## YOUR CODE HERE
# region <- getData(country = 'DNK', level = 1)
region <- readRDS("../data/gadm36_DNK_1_sp.rds")

# select Midtjylland
Midt <- region %>% as("sf") %>% filter(VARNAME_1 == "Central Jutland")

# intersect with municipalities
no_cars_pct_midt <- no_cars_pct_mun %>% st_transform(4326) %>% st_intersection(Midt$geometry)

# quick look that we got the right area
mapview(no_cars_pct_midt, zcol = "nc2021pct")

Gosh, that was a ton of wrangling! But we are nearly there! :)

Spatial overlay with sf

Spatial overlay is a very common operation when working with spatial data. It can be used to determine which features in one spatial layer overlap with another spatial layer, or extract data from a layer based on geographic information. In R, spatial overlay can be integrated directly into tidyverse-style data analysis pipelines using functions in sf. In our example, we want to determine the areas in Midtjylland with the greatest proportion of households without access to a car that also are beyond a 20 minute walk/cycling route from a hospital.

To do this, follow these steps:

  • doublecheck that the coordinate reference system of the no_car_pct_midt dataset is 4326, the same CRS used by the isochrones;
  • extract only those municipalities with a percentage of households without cars of 25 percent or above;
  • use the st_difference() function to “cut out” areas from those municipalities that overlap the 20-minute cycling isochrones. (Union the isochrones first for neater output).
  • Once you complete this operation, visualize the result in your Mapbox map.
## YOUR CODE HERE

Questions:

3. Are there any areas that are located beyond a 20-minute bike-ride of a hospital that have proportionally lower access to cars? Where are they and what is their total area?

4. How would you ensure residents in these areas have reasonable access to hospitals?

target25pct <- st_union(target_areas)  # gets you area in sq kms
mapview(target25pct)
st_area(target25pct)/1e+06
7968.367 [m^2]